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D , We calculate various thermodynamic quantities of vortex liquids in a layered superconductor 

' by using the nonperturbative parquet approximation method, which was previously used to study 

the effect of thermal fluctuations in two-dimensional vortex systems. We find there is a first-order 
transition between two vortex liquid phases which differ in the magnitude of their correlation lengths. 
As the coupling between the layers increases, the first-order transition line ends at a critical point. 
We discuss the possible relation between this critical end-point and the disappearance of the first- 
C ■ order transition which is observed in experiments on high temperature superconductors at low 

O ' magnetic fields. 

■ 
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1 , I. INTRODUCTION 

Q ' Thermal fluctuations and quenched disorder are two important factors responsible for the rich phase diagram 
O I exhibited by high temperature superconductors in a magnetic field . Recent experiments [^jsj on high temperature 
superconductors show a first-order transition line well below the upper critical field Hc2{T) in the field-temperature 
^-H \ (H-T) phase diagram. When the strength of disorder is weak, this line is usually interpreted as a melting line of a 
J> • vortex lattice into a vortex liquid. The first-order transition disappears at both high and low magnetic fields 

\ In the vortex lattice melting picture the loss of first-order behavior at high fields is usually attributed to the effect of 

■ disorder, which is supposed somehow to produce a multi-critical point where the first-order transition changes into a 
continuous one. 

In a recent numerical simulation Q of a vortex system in a layered superconductor, a first-order transition line was 
^) \ also obtained which disappeared at a critical end-point at low magnetic fields. According to the simulation results, 

■ in contrast to the vortex lattice melting picture, there is only one phase below and above the transition, namely the 
vortex liquid but with different correlation lengths on either side of the first-order transition line. However, in other 
simulations using periodic boundary conditions within a layer, an apparent first-order transition between vortex 

^ liquid and lattice was obtained. 

In this paper we apply an analytical approach to a layered superconductor in a magnetic field perpendicular to the 
P5 layers in an attempt to elucidate the nature of phase transitions in the system. We use the parquet approximation 
O which has been successfully applied to the two-dimensional vortex system. It is a nonperturbative analytic 

method free from any finite size or boundary effects in the direction perpendicular to the magnetic field. The parquet 
approximation deals with the renormalized four-point function of the vortex system which is obtained by summing an 
infinite subset of Feynman diagrams, the so-called parquet diagrams. In two dimensions this method is able to capture 



X 



the growing crystalline order developing in the vortex liquid as the temperature is lowered |10|. As explained in the 
I discussion below, the introduction of an additional dimension to the parquet calculations imposes severe numerical 
difficulties because of the large number of variables specifying the renormalized vertex function. This problem has 
restricted us in this paper to considering a small system which consists of just four layers satisfying a periodic 
boundary condition along the field direction. Each layer is, however, treated as an infinite plane so we can are in 
effect descibing an infinite number of vortices. Despite the unphysically small size of the system we consider here, we 
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find that the efi^ect of the inter-layer couplng makes a significant difference compared with the two-dimensional case, 
where no finite-temperature phase transition occurs in the parquet calculation. In the layered system, we find that 
the thermodynamic quantities describing the vortex liquid show abrupt changes when the field and the temperature 
are varied, and these sharp changes are interpreted as a first-order phase transition within the vortex liquid phase. 
We find that the first-order transition line ends at a critical end-point. The two phases we obtain above and below the 
transition are both liquid phases with different correlation lengths, consistent with the numerical simulation results 

In the next section we present our model for the vortex liquid in a layered superconductor. We then set up the 
parquet equations for the layered system. In the final section we present the numerical solutions of the parquet 
equation for the four-layer system which are interpreted in terms of a first-order transition and the termination of 
the first-order transition line at a critical end-point. We conclude with a discussion on the possible implication of our 
results for the H-T phase diagram of a layered superconductor. 

II. MODEL 

Our starting point is the Lawrence-Doniach model for a layered superconductor in a magnetic field perpendicular 
to the layers. With the order parameter in the n-th layer denoted by ipn, the free energy is given by 
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where do is the layer thickness, d the layer spacing and a, /3, mat, and rUc phenomenological parameters. We denote 
by T = h'^/2mcd'^ — {^c/d)'^ the dimensionless ratio between the coherence length <^c perpendicular to the layers and 
the layer spacing. We take B = V x A as constant and uniform. 

We use the lowest Landau level (LLL) approximation which is believed to be valid over a large portion of the vortex 
liquid regime. We expand the order parameter '0„(r) in terms of the eigenstates of \{~ifiW — ^A)-i/;„p and keep 
only the lowest eigenvalue state. In the symmetric gauge, where A = ?/, x, 0), the LLL wavefunction is given by 
^/jlf"^{r) = exp(— ^^|z|^/4)(/)„(z) where fj,"^ = e*B/hc and 0„(z) is an arbitrary analytic function oi z = x + iy. In the 
LLL approximation, the free energy becomes 

=^do f dz*dz[aHe-^''\'-\'/^\Mz) 

+ re-^'l-l'/2|,/,„(z)-0„+i(z)|2], (2) 

where an = a + e* BTi/2cmab changes sign crossing the upper critical field line Hc2{T). Physical proper- 
ties of the vortex system in the layered superconductor are then determined from the partition function Z = 
/ n„ 2?<^n2?0; exp(-f^[0, 0*]/fcBr). 

For a two-dimensional system, that is for a single layer, the temperature and field dependence are all contained 
within the single dimensionless parameter a2T = (27rdo//3M^^s^)^^^Q!-H- The inter-layer coupling strength is described 
by the dimensionless parameter tt = (27rdo//3/i^fcsr)^/^T. For the layered system it is very convenient to use the 
dimensionless field-temperature parameter ax employed in studies of the continuous anisotropic three dimensional 
GL model: ar — (S-kHc/ (3'e* BkBT)'^/'^{h^ /2mc)^^'^aH where f3' — {d/do)P is the coefficient of the quartic term for 
the three dimensional GL model. The two-dimensional limit corresponds to 0, while in the limit — > oo for 

constant ay, the system behaves as a continuous three dimensional model. Note that ar = 2'^^^T^^^a2T- In terms of 
t = T/T,o aTidh = H/HMO), 
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where Tco is the zero- field transition temperature and Hc2 (0) is the straight line extrapolation of the Hc2 (T) line near 
Tco to zero temperature. 

We shall calculate various correlation functions using the parquet approximation. This is a nonperturbative analytic 
approximation and requires no boundary conditions to be imposed in the direction perpendicular to the magnetic field 
within each layer. Our system consists of a stack of N layers, on which we impose the^periodic boundary condition 
in the direction along the field 4>n+N{z) = 4>n{,z). We introduce the Fourier transform of (/>„ via 

= E (4) 

qelst.B.Z. 

where q — 271] /Nd and we use N integer values of j in the sum such that the wavevector q e [— Tr/d, Tr/d) belongs to 
the first Brillouin zone. 

One can develop the standard perturbation theory from the given partition function. The bare propagators are 
given by 

(0;,(z'*)0,(z))o = <eN5,.,,{p^/2^)e^^''"^/^G^{q), 

where 

^'^'^^ ^ (^) aH + 2T{l-cos{qd)y 

Since the magnetic length /i^^ is the only length scale perpendicular to the field direction which appears in the 
propagator pO[ |, the fully renormalized propagator can also be written as (|^) with the renormalized G{q) replacing 
the bare function Go{q). The renormalized G{q) is determined self-consistently in the parquet approximation. 
^The niain quantity one studies in the parquet approximation is the renormalized connected four-point function, 
('/'gi (-2^i)^^q2 ('^2)'^<?3 (^3)'/'g4 (•^4))c which can be written in terms of the renormalized vertex function r(qi, 92, ^a; k) = 
r(9i, 92, 93, (Zi +92 — Qa; k) Here the wavevector k lies in the two dimensional space perpendicular to the magnetic 
field and qi is a wavevector along the field direction. Note that to the lowest order r(9i, 92, 93; k) — TbO^-) independent 
of the wavevectors along the field direction. 

In order to make a resummation over all parquet diagrams, we note that the contributions to F can be decomposed 
into a totally irreducible part and a reducible part which in turn can be written as the sum of three parts = 
1,2,3) representing the contributions from three different channels. (A detailed discussion on the diagrammatic 
decomposition can be found in Ref. ||l0|] ). The parquet approximation we employ here corresponds to neglecting in 
the totally irreducible vertex all the higher order (O (/?'*)) diagrams except the bare vertex function FB(k) so 

3 

r(9i,92,93;k) = FB(k) +^nj(9i, 92,93; k). (6) 

1=1 

Each reducible vertex Hi is composed of the irreducible vertex where 

A,;(9i,92,93;k) = rs(k) + ^ nj(9i, 92, 93; k) (7) 

and the renormalized F via the following Bethe-Salpeter equations: 

ni(9i,92,93;k) = -^^G{p)G{qi +92 -p) Ai(9i, 92,_p) o F(p, 91 + 92 -^,93) (k), 
p 

n2(9i,92,93;k) = -—^G{p)G{p-qi + 93)A2(9i,P - 9i + 93, 93; k)F(p, 92,_p + 93 -9i;k), 
p 

2 ^ ^ r -1 

n3(9i,92,93;k) = -— ^G(p)G(p + 92 - 93) A3(9i,p + 92 - 93, p) * F(p, 92, 93) (k), (8) 
p 

where the operation o between two arbitrary functions /(k) and ^(k) is defined by 

(/o5)(k) = -y ^/(k-k').9(k') 

X cos{{k^k'y - kyk'J/fi'^), 
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and f * g is just the convolution which is the same as the above expression without the cosine term. In (||) we have 
used the dimensionless form of the propagator function G{q) = {d{)l3y? /2'KkBTY/'^G{q). The above equations, (|^), 
(0) and ^ form a closed set for V when the renormalized propagator G{q) is known. In the parquet approximation 
G{q) is determined self-consistently by the Dyson equation: 



G-\q)~G^\q)^-Y^G{q') 



^ G{q')G{q")G{q + q'-q")^ J ^r{q, q' , q" ;k 

q\q" 



(9) 



Using the solutions to the above equations one can calculate several interesting physical quantities. The structure 
factor, which is the measure of correlation between vortices in the vortex liquid, is calculated from 

X„-„'(r-r') = (|*„(r)P|vI/„,(r')P) 

-(|vl/„(r)p)(vl/„,(r')p). 

The structure factor A„j(k) used in this paper is then defined by 

A™(k) ^ {^y'''' J d^Re*-^x™(R). (10) 
Using (^, one can write the generalized Abrikosov ratio Pa as 

{\^n{r)n _^-N-'EqGiq)G,\q) 



Ha 



(l*n(r)P)2 



(11) 



^-^Y.qG{q) 

The Josephson coupling parameter 77 measures the degree of independence between adjacent layers and is given by 

|M/„(r) - «'n+i(r)n _ iV-i 2(1 - cos(gd))G(g) 



(l*«(r)|2) iV-iE^G(g) 
We can also put the above equations in more convenient form as follows. If we introduce 

^ ^ ^ -11/2 
7(91, 92, 93; k) = G{qi)G{q2)G{q3)G{qi + q2 - qs) r(qi , ^2, 93; k). 



(12) 



(13) 



and similary Ai((7i, (72, 93! k) and tt^ (gi , (72 , 93 ; k) from Ai(gi, (72, 93; k) and n,;((7i, 1721 93! k) respectively, and denote 
the Fourier transform of 7, n and A by 7({ni};k), 7ri({ni};k) and Ai({ni};k) with the shorthand notation {rii} = 
(711,77.2,713,714), then (@) is simplified to 



7ri({ni};k) = - ^ ^Xi{ni,n2,n' ,n") o j{n' ,n" ,n3,ni) (k) 
7r2({77i};k) = -2 ^ A2(77,i, n', 7^3, t^"; k)7(77,", 712, 7x', 714; k) 

n' ,n" 

({?7i};k) = ~2 ^ |^A3(77,i, 77,', 7i", 7^4) * 7(71", 77,2, 77,3,71') (k). 



The remaining parquet equations can easily be derived as follows: 

Aj({?7,i};k) = 7({77j};k) - 7rj({77j};k) 

for j — 1, 2, 3 and 

3 

7({nJ; k) = fB{{n,})TB{-k) + ^,({r7.}; k). 



(14) 



(15) 



(16) 
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where 



fB{{ni]) ^Y.^{1 - n^)g{l - n2) 
I 

xg{n3-l)gin4-l) (17) 
with g{n) being the Fourier transform of G^/^{q): Q{n) — (1/-/V) '^^ey:.j>{iqnd)G^/'^{q). 

III. RESULTS 

We solved (numericaUy!) the parquet equations, (|l^-(|l6|) for 7({ni};k) with Eq. (||) for the propagator G{q). 
One could use equivalently the set of equations, (H)-^ for r({(ji}; k). It is just a matter of convenience. In the 
present analysis the former was used. In solving we start from some initial functions, Ai, 7 and G{q) and update 
these functions iteratively using ([l^)-([l^) and (^). As mentioned above, the main numerical difficulty compared to the 
two-dimensional case is that 7 and A.; now involve three extra indices in addition to the two-dimensional momentum k. 
A large amount of computer memory and CPU time is required as the number of layers becomes large. In this paper 
we only consider a small system where the number of layers is four and try to see what effect the coupling between 
the layers has on the two-dimensional parquet results. At high temperatures we find that the iteration converges very 
quickly, but as the temperature is lowered the convergence gets slower, and furthermore we have to use the results 
obtained at a nearby temperature as the initial values for A^ and 7 in order to get rapid convergence. 

From the solutions F and G{q) to the parquet equations for given a2T and tt, we calculated the thermodynamic 
quantities introduced in the previous section. The three-dimensional temperature parameter aT can be obtained from 
the two parameters. Figures 1 and 2 show the Abrikosov ratio (3a, defined in Eq. ( [Til ) and the Josephson coupling 
parameter rj in Eq. ( p^ ) as functions of ut for various values of the interlayer coupling tt- We find that for some 
values of t^, these thermodynamic quantities show abrupt changes as q^t is varied. But this behavior disappears 
for Tt > 0.11 when now j3A and r/ just decrease smoothly without showing any sudden drops as the temperature is 
lowered. These features can be explained by the existence of a first-order transition in the vortex liquid and of a 
critical end-point: the first-order nature of the transition gets weaker as the inter-layer coupling strength tt increases, 
eventually disappearing at the critical end-point at ~ 0.11. 
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FIG. 1. The Abrikosov ratio /?a as a function of temperature parameter ut for different inter- layer coupling strengths 
tt = 0.1 (filled squares), 0.105 (open squares), 0.11 (open triangles) and 0.12 (filled triangles). The dotted lines are guides to 
the eye. 




FIG. 2. The Josephson coupling parameter r/ as a function of temperature parameter qt for different inter-layer coupling 
strengths r. The symbols are the same as in Fig. 1. The dotted lines are guides to the eye. 



The sudden changes in these quantities occur in a very narrow range of the temperature parameter ay. We believe 
it is not a crossover which just happens to resemble a first order transition, since we find instabilities in obtaining the 
solutions to the parquet equations within this narrow range of (Xt- This is illustrated in Fig. 3. It shows a typical 
example of how thermodynamic quantities such as (^a or ry evolve as we itcratively solve the parquet equations in 
the vicinity of the transition. Well above the transition, (3a converges rather quickly as can be seen in Fig. 3. This 
is also the case well below the transition. However, as the transition is approached from above the convergence gets 
slower, then we reach a very narrow region of where the iteration appears not to be converging. We find this kind 
of instability for the cases where tt < 0.11, while no instability is seen for tt > 0.11. This is in fact how we locate 
the critical end-point in our model. Since the iteration method works well only if we use initial functions which are 
close to the actual sohitions, we conclude that this instability signals a sharp change emerging in the system. Strictly 
speaking, there should exist a solution for any value of o.t even if there is a first order transition. Therefore the 
instability does not imply that there is no convergent solution. It just tells us that it is very hard to get a solution 
using the iteration method employed here. (It is interesting to note that the solution below the transition was in 



6 



fact found after a very long iteration). For this reason, although we find that the size of the first-order transition- like 
jumps in thermodynamic quantities such as Pa decreases as the critical end-point is approached, it is quite hard to 
determine accurately the size of these jumps. 



1.55 
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FIG. 3. The evolution of the Abrikosov ratio /3a in the iteration of the parquet equations. The inter-layer coupling is fixed 
at tt ~ 1.075. The three dimensional temperature parameter qt is given by -3.594 (solid), -3.654 (dashed), -3.684 (dotted), 
and -3.690 (dot dashed). 

The first order transitions we obtain here are between two vortex liquid phases with the phase below the transition 
being more correlated. This can be seen in Fig. 4 where the structure factor An{K) defined in Eq. (|l^) is shown at 
temperatures above and below the transition. In both cases, there is no crystalline long-range order present in the 
system. We note that the first peak at ~ 2.69, which is the value of the first reciprocal lattice vector of a triangular 
lattice, represents the crystalline order developing in the plane perpendicular to the magnetic field. Comparing the 
structure factors above and below the transition, the difference is most significant for the case of n = 2. This tells us 
that the state below the transition is a vortex liquid with larger correlation along the field direction. Therefore the 
transition we find is more like a decoupling transition occuring in the vortex liquid. This can also be seen in the drop 
in r/ (Fig. 2) as the temperature is lowered, since rj is largest when neighboring layers are more or less independent. 



(a) (b) 




K 



FIG. 4. The structure factor An{K) at tt = 0.105 and four different values of ar- The dimensionless two-dimensional 
wavevector K = k/jj.. (a) n=0,(b) n=l, and (c) n=2. For all cases, aT= -3.715, -3.709, -3.685 and -3.566 reading from top to 
bottom. 
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Figure 5 shows a collection of the first-order transition points in the ut-tt plane and the location of the critical 
end-point. The transition line is almost independent of ax- Although we were only able to study the region close to 
the critical end-point, the nearly aT-independent transition line is consistent with the result obtained in Ref. ||]. We 
found it very hard to extend this line to the small- region. Since the size of the jump gets bigger as we move away 
from the critical end-point to the small- tt region, a solution below the transition is going to be very different from 
the one above, which makes the parquet equations hard to solve by iteration. According to Eq. (||) the transition line 
which ends at large tt corresponds to a critical end-point at low magnetic field. In the H-T space, the shape of the 
phase boundary looks like the one found in Ref. Q, which is consistent with the experimental results on YBCO-type 
superconductors. However, we do not expect that the actual position of the transition line in our work can be directly 
compared with experiment since numerical difficulties have limited us to studying only a small number of layers. The 
transition temperatures {aT ~ —3.7) in our model are higher than those found in Ref. |^] (ax ^ —7.0), and therefore 
are also higher than the experimental values. 
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FIG. 5. The first order transition line and the location of the critical end-point (cross) in the ot-tt space. The inset shows 
how the size of the jump in /3a at the transition decreases as the critical end-point is approached. 

To summarize we have applied a nonperturbative analytic method to the vortex liquid system in a layered supercon- 
ductor. The inter-layer coupling produces a first-order transition line which ends at a critical end-point at low fields, 
whereas for a purely two-dimensional system there are no transitions of any kind within the parquet approximation. 
The results are consistent with the first-order transition being a decoupling transition between two vortex liquid 
phases. Clearly in order to extend this method to other situations eg. lower temperatures, the effects of disorder and 
above all, more layers, we have to devise a way to speed up the convergence of the iterative solution to the parquet 
equations. One possible method might be to use a combination of solutions obtained in previous steps as the next 
stage of the iteration . 
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